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We explore the properties of the bipolaron in a ID Holstein-Hubbard model with dynamical quantum 
phonons. Using a recently developed variational method combined with analytical strong coupling 
calculations, we compute correlation functions, effective mass, bipolaron isotope effect and the phase 
diagram. The two site bipolaron has a significantly reduced mass and isotope effect compared to the 
on-site bipolaron, and is bound in the strong coupling regime up to twice the Ifubbard U naively 
expected. The model can be described in this regime as an effective t-J-V model with nearest 
neighbor repulsion. These are the most accurate bipolaron calculations to date. 



PACS: 74.20. Mn, 7f.38.4-i, 74.25. Kc 

While there is a generally accepted belief that in high 
Tc superconductors a dominantly electronic interaction 
is responsible for the unusually high transition temper- 
atures, the interplay between the electron-phonon inter- 
action and the strong electron-electron interaction nev- 
ertheless plays a significant role in determining the phys- 
ical properties of these strongly correlated systems . 
Although the study of lattice effects in strongly corre- 
lated materials is steadily growing 1^-0 , the understand- 
ing of the influence of the Hubbard interaction U on 
bipolaron formation is still incomplete. In particular, it 
is known that in the strong coupling regime bipolarons 
have an extremely large effective mass ||,^ , which repre- 
sents one of the main objections |^Jl^ against the theory 
of small bipolaron superconductivity |^]. Recent calcu- 
lations in the adiabatic (static phonon) limit show that 
a first-order phase transition exists between the on-site 
(SO) and lighter neighboring-site (SI) bipolaron ||l^]. The 
properties of the SI bipolaron were also investigated by 
variational and exact diagonalization methods |l^,|lj] . 

In this letter we present accurate numerical solutions 
of the Holstein-Hubbard model bipolaron on a ID infinite 
lattice. Our results are exact to within the line-widths 
on the figures in the intermediate and near the strong 
coupling regimes. The results are compared to analytical 
calculations in the strong coupling regime. 

We consider the Holstein-Hubbard Hamiltonian flSl 
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where ^ creates an electron of spin s and creates a 
phonon on site j. The last term in Eq. (|l|) represents the 
on-site Coulomb repulsion. We consider the case where 
two electrons with opposite spins couple to dispersionless 
optical phonons. 

Basis states for the many-body Hilbert space can be 
written \M) = ^2; ■ ■ . ,n„,nm-i-i, •■■,), where the up 



and down electrons are on sites jl and j2, and there are 
rim phonons on site m. A variational subspace is con- 
structed beginning with an initial state where both elec- 
trons are on the same site with no phonons, and operating 
repeatedly (A'^^.-times) with the off-diagonal pieces {t and 
g) of the Hamiltonian, Eq. (0). All translations of these 
states are included on an infinite lattice. This method 
was used previously for the polaron (one electron) |]l6| , |l7| . 
The method is very efficient in the intermediate coupling 
regime, where it provides results that are variational in 
the thermodynamic limit and bipolaron energies that are 
accurate to 7 digits for the case Nh — 18 and size of the 
Hilbert space N = 2.2 x 10^ phonon and down electron 
configurations for a given up electron position. 

Before presenting the numerical calculations, we show 
that many interesting properties of the bipolaron can be 
found in second-order strong coupling perturbation the- 
ory. Following Lang and Firsov we use the canonical 
transformation H — e^He~^, where S = 
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The transformed Hamiltonian takes the form 



H = Ho + T 
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where 
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The first term in Hq is the en- 



ergy of the phonon excitations, and the second is the en- 
ergy gained by displacing the oscillator in the force of the 
electron. The exponential factors in the T hopping term 
arise because the Lang-Firsov transformation redefines 
the origin of the harmonic oscillator when the number of 
electrons on the site changes. In the limit g ~> and 
uj ~* 00 with ujg^ constant, the phonon interaction is in- 
stantaneous and the Holstein-Hubbard model maps onto 
a Hubbard model with an effective Hubbard interaction 
U — U ~ "iuig^. In ID a bound state exists for U < 0. 

In the strong coupling or antiadiabatic limit, T in 
Eq. (0) is considered a perturbation. It represents the 



1 



hopping of electrons, including possible creation and de- 
struction of phonon excitations. The SO state (f)o = 
cj|cjj0) has the lowest energy to zeroth order in T when 

U < 0. In perturbation theory to second order, the en- 
ergy of the SO bipolaron is 
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The effective mass m^g = d"^ Ebi{k) / dk"^ can be obtained 
from Eq. (p|) by calculating the infinite sums [p| 
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where r(a;) and V{a,x) are the gamma and incom- 
plete gamma functions respectively. Taking the large 
g limit in Eq. (P, one finds (see also Ref. |^) m'^Q ~ 
4-/??^^ exp(— 45^)/(w5) for [/ = 0. At large g, mso is 
roughly a factor exp(3(7^) larger than the polaron mass, 
m~o - 2texp{-g^). 

One would naively expect that within the strong cou- 
pling approximation, a bipolaron unbinds when U > 0. 
This is false: a bound bipolaron may exist even for U > 0. 
In this regime, a set of degenerate states (/>, 
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for i ^ 0, written in a translationally invariant form, rep- 
resents states with minimum energy of Hg. The energy of 
a SI bipolaron is obtained by solving the secular equation 



0, where matrix elements = {(l)i\T\(l)j) 



are calculated to second order in T. 

The main source of binding arises because the second 
order matrix element Ti_i for two neighboring electrons 
Iti) to exchange sites is exponentially larger than all 
other off-diagonal matrix elements. While in the limit 
g oo, Ti_i cx t'^/U, the second largest (first order) 
matrix elements T^^^+i c>c iexp(— g^) for i ^ {0,-1}. 
The exponential suppression arises from phonon rear- 
rangement whenever the initial and final charge distri- 
butions differ. In the singlet configuration, 0f"° = 
(01 -I- (/)_i)/\/2, the diagonal correction to the energy is 
given by Tn -I- Ti_i. There is also a contribution to Tn 
that resembles a retardation effect, in which one elec- 
tron hops creating one or more phonon quanta on its 
initial site, and the second electron follows absorbing the 
phonons. This process, however, decreases exponentially 
with g as exp{— g'^) / {ujg'^) , and is not strong enough to 
bind two polarons in the triplet configuration. 

To obtain the SI bipolaron stability criterion, we note 
that the secular equation for fixed center of mass mo- 
mentum k can be mapped onto a tight-binding model 
for a linear chain. Since all the off-diagonal matrix ele- 
ments except Ti_i scale as exp(— g^), the condition for a 
bound state is given by Tn -I- Ti-i < Tu. Keeping only 



terms of order t^/{u}g'^) and setting /c = 0, the matrix 
elements can be expressed as Tu = + 7],, Ti_i = Tb, 
and Tii = 2Ta\ i > 1, with 
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where the final terms refer to the large g limit. In the 
large g hmit, binding occurs for U < Aujg"^. 

The effective SI bipolaron mass is computed by ap- 
proximating the SI bipolaron wavefunction with only 



(omitting the exponential tail), 
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There are three distinct processes contributing to msi- 
an SI pair can move by one lattice site through an in- 
termediate doubly occupied state with n phonons (terms 
that contain U in Eq. (0)), or through an intermediate 
state with only phonon degrees of freedom (terms with- 
out U). 

The bipolaron isotope effect is a measure of how 
the bipolaron mass varies with the ion mass M, 
a = c?lnmbi/c?ln M . This can also be written a — 
— ^ dlnnihi/d In LU, where the derivative is taken with 
Lug'^ held constant. The bipolaron isotope effect a is 
equal to the superconductivity isotope effect asc = 
— d In Tc/d In M only in the low density limit, when super- 
conductivity occurs as a weakly interacting gas of bipo- 
larons Bose condenses. Then in mean field theory (ignor- 
ing fluctuations that are important in low dimensions), 
the transition temperature Tc at fixed density is propor- 
tional to TO(~^, and asc = ct- The bipolaron isotope effect 
is expected to change substantially between the SO and SI 
regimes. Indeed, using Eqs. (^ ^ in the strong coupling 
limit, we obtain aso ^ 2g^ — j for U — and aso de- 
creasing with increasing U. In the SI regime, asi ''^ 3^/2 
and is only weakly dependent on U. The bipolaron iso- 
tope effect has also been calculated by Alexandrov [p^ ; 
the superconductivity isotope effect has been measured 
by many groups pO|| . 

We now present numerical variational results, using 
units where hopping t — 1. Fig. (0) shows the ground 
state electron-electron density correlation function C{i — 
j) = (V^hi"-jlV'o), where = -|- rii^. At 5 = I, 
the bipolaron widens with increasing U and transforms 
into two unbound polarons (which can only move a finite 
distance apart in the variational space). The value U — 
1.5 is below the transition to the unbound state at Uc — 
2.17, calculated by comparing the polaron and bipolaron 
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energies. At U = 1.5, C{i) falls off exponentially, while 
for U > Uc the typical distance between electrons is the 
order of the maximum allowed separation Nh A 
state of separated polarons is clearly seen for U — 20. 
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FIG. 1. Electron-electron correlation function C{i — j) cal- 
culated at cj = 1, a) (/ = 1 and b) <; = 2 for different values 
of U, with Nh = 18. The two ordinate axes have a different 
range. Insets show results for U = 0. All curves are normal- 
ized, = 1. 

Two distinct regimes are seen at g — 2 within the 
bipolaronic region. At U — 7 < Uq = lug^ = 8, the 
correlation function represents the SO bipolaron, while 
aX U = 9 > C/o we find the largest probability for two 
electrons to be on neighboring sites, which is character- 
istic of the SI bipolaron. In contrast to previous calcu- 
lations where phonons were treated classically (l2|, we 
find a crossover rather than a phase transition between 
the two regimes. The precision of presented correlation 
functions in the bipolaron regime is within the size of the 
plot symbols in the thermodynamic limit. 

Figure (||a) plots the bipolaron mass ratio Rm = 
mu/2mpo vs. U for different values of lo and g. In all 
cases presented in Fig. (^), Rm approaches 1 as [/ ap- 
proaches U = ?7c in agreement with a state of two free 
polarons. At fixed lo = 1 the bipolaron mass ratio in- 
creases by several orders of magnitude with increasing g 
at [/ = 0. The increase can be understood within strong 
coupling theory. Increasing U sharply decreases Rm in 
the SO regime. Note that the scale in Fig. ^ is logarith- 
mic. Near the strong coupling regime {g = 2) and for 
small J7, good agreement is found between the numeri- 
cal and the strong coupling result obtained from Eq. (|^). 
The difference between these results increases as U ap- 
proaches Uq = 8, where the perturbation theory based on 
the SO bipolaron breaks down. In the SI regime U > Uq, 
Rm is small, as predicted by the strong coupling result. 
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FIG. 2. a) The mass ratio Rm = mbi/2mpo vs. U and 
b) the bipolaron isotope effect a vs. U. Numerical results 
are for Nh = 18. Results for Rm at ui — 0.5 are obtained 
by extrapolating A^'^ oo. Precision in all curves is within 
the line- width in the thermodynamic limit, except for a with 
Lj — 0.5, where the error is estimated to be ±5%. The thin 
line in a) is the strong coupling expansion result for uj = 1, 
g = 2. Polaron masses in units of the noninteracting elec- 
tron mass are nipo = 1.35, 1.76, 10.4, 3.06 from top to bottom 
respectively. 

The bipolaron isotope effect, shown in Fig. (^b), is 
large in the strong coupling (w — l,g — 2) and small 
U regime, where its value is somewhat below the large 
g strong couphng prediction aso ^ 2g^ — i = 7.75. 
With increasing U, a decreases and in the SI regime 
approaches asi = 5^/2 = 2. A kink is observed in the 
crossover regime. 

We conclude with the phase diagram Uc{g) shown in 
Fig. (H) at fixed to = 1. Numerical results, shown as 
circles, indicate the phase boundary between two dis- 
sociated polarons each having energy Epo and a bipo- 
laron bound state with energy Ei,i. In the inset of 
Fig. (j^) we show the bipolaron binding energy defined 
as A = Ehi — 2Epo. The phase diagram is obtained 
from A = 0. The dashed line, given by Uq = 2ojg^, is 
a reasonable estimate for the phase boundary at small 
g. At large g the dashed line roughly represents the 
crossover between a massive SO and lighter SI bipolaron. 
The dot-dashed line is the phase boundary between SI 
and the unbound polaron phase, as obtained by degen- 
erate strong coupling perturbation theory. Numerical re- 
sults approach this line at larger g. The dot-dashed line 
asymptotically approaches Uc — Aiog"^. 

In the SI regime, the strong coupling expansion can 
be used to rewrite the Holstein-Hubbard Hamiltonian 
Eqs. (|l|,§) as an effective t-J-V model that applies to an 
arbitrary number of particles. 
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Heff = -h ^ c|,sCi+i,s + h.c. + ^ JSiS,j+i + Vnin.,+i, 

i,s i 

where h = i exp(-52), J = -2Tt, > 0,V = -Ta+Tb/2 > 
and Ci_s — Ci^si^ — ni^-s) p2| . The 1^ term is a repulsion 
between nearest neighbors. For simplicity and in keeping 
with the approximations used to derive the standard t—J 
model, we have omitted next-nearest neighbor hopping 
terms that are of order exp{— g^) / (ujg^) , and a constant 
term. While the standard t — J model can be derived 
from the Hubbard model only for J < t, the parameters 
of Heff are quite different, with J typically much larger 
than ti. In the static limit, an SI bipolaron is stable in 
the interval J/4 < V < 3 J/4. In the dilute limit, only 
singlet bipolarons exist with binding energy (in the static 

limit) ^5=0 = -3 J/4 + V At^/U + t^/{ujg^), while 

triplet fluctuations are almost completely frozen out due 
to the high energy i?s=i — J/'i + V- Such a state is 
very similar to a negative-C/ Hubbard model, except that 
singlets occupy two sites, like in the RVB model. 



the Holstein-Hubbard model, where antiferromagnetism 
is a consequence of the original Hubbard interaction U, 
while the attractive interaction is mediated by phonons. 
Taking into account the similarity between a system of 
SI bipolarons and the negative-U Hubbard model, one 
should not rule out the possibility that SI bipolarons 
form a superconducting state of either s— wave or d— wave 
symmetry in two spatial dimensions. In such a state, 
strong electron-electron interactions and electron-phonon 
coupling should be treated on an equal footing. 
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FIG. 3. Phase diagram and binding energy A in units of t 
(inset) calculated at = 1. Numerical results are circles. For 
greater accuracy, results near the weak and strong coupling 
regime were obtained by extrapolating A*';, oo. 

In conclusion, we have demonstrated that near the 
strong coupling limit a mobile SI bipolaron exists with 
an effective mass of the order of a polaron mass, and 
an isotope effect in the range < a < The 
wavefunction of the SI bipolaron is a spin singlet with 
extended s— wave spatial symmetry. Taking into ac- 
count the asymptotic stability criterion Uc = Aojg'^, it is 
clear that a triplet SI bipolaron that corresponds to the 
U —^ oo solution is not bound. In the static limit, it can 
be shown that bound states of three or more polarons are 
not stable in the SI regime, thus ruling out phase sep- 
aration. An effective t-J-V model captures the physics 
of many SI bipolarons in the strong coupling regime of 
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